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Abstract 

The Macaulay2 package GraphicalModels contains algorithms for the algebraic study of graphical 
models associated to undirected, directed and mixed graphs, and associated collections of conditional 
independence statements. Among the algorithms implemented are procedures for computing the vanish- 
ing ideal of graphical models, for generating conditional independence ideals of families of independence 
statements associated to graphs, and for checking for identifiable parameters in Gaussian mixed graph 
models. These procedures can be used to study fundamental problems about graphical models. 

1 Graphical models 

A graphical model is a statistical model associated to a graph, where the nodes of the graph represent random 
variables and the edges of the graph encode relationships between the random variables. Graphical models 
are an important class of statistical models used in many applications (see standard textbooks [5, 8 ) because 
of their ability to model complex interactions between many random variables, by specifying interactions 
using only local information about connectivity between the vertices in a graph. 

There are two natural ways to specify a graphical model, through either conditional independence state- 
ments specified by the graph or via a parametric representation (often called a "factorization"). Every 
distribution which factors according to the graph satisfies the conditional independence statements implied 
by the graph. This leads to the question: Which distributions satisfy the conditional independence state- 
ments implied by the graph, but do not factor? 

Once we specify the types of random variables under consideration (e.g. discrete random variables 
or Gaussian random variables) it is possible to address the questions in the preceding paragraph using 
(computational) algebraic geometry. Indeed, in these cases, the set of all probability distributions satisfying 
a family of conditional independence constraints is a semialgebraic set. For discrete random variables, that 
semialgebraic set is a subset of the probability simplex, and can be represented by a certain homogeneous 
ideal generated by quadrics. For Gaussian random variables, this set of distributions corresponds to a 
semialgebraic subset of the cone of positive definite matrices. Similarly, the parametrized family of probability 
distributions also is a semialgebraic set (of the probability simplex for discrete random variables, and of the 
cone of positive definite matrices for Gaussian random variables). This algebraic perspective has been studied 
by many authors [31H1IS], and the book pQ provides many details. 

The Macaulay2 package GraphicalModels allows the user to compute the ideals of conditional indepen- 
dence statements for any collection of statements for discrete or Gaussian random variables. It can also 
compute the vanishing ideal of a graphical model in these cases. A number of auxiliary functions are useful 
for doing further analyses of graphical models. 

To be explicit, consider the directed acyclic graph G with five vertices {a, 6, c, d, e} and edge set {a — > 
d, b — > g?, c — > d, c — ¥ e, d — ¥ e}. The following commands compute the associated conditional independence 
ideal for the set of global Markov statements, CJ g i bai(G) i an d the vanishing ideal Iq of the Gaussian graphical 
model on G. 
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Macaulay2, version 1.4 

11 : loadPackage "GraphicalModels" ; 

12 : G = digraph{{a ) d},{b ) d},{c,{d,e}},{d,e»; 

13 : R = gaussianRing G 

03 : PolynomialRing 

14 : vars R 

04 = I s_(a,a) s_(a,b) s_(a,c) s_(a,d) s_(a,e) s_(b,b) s_(b,c) s_(b,d) s_(b,e) s_(c,c) 

s_(c,d) s_(c,e) s_(d,d) s_(d,e) s_(e,e) I 

14 : I = conditionalIndependenceIdeal(R,globalMarkov(G) ) ; 

15 : J = gaussianVanishingldeal(R) ; 

16 : flatten degrees J 

06 = {1, 1, 1, 2, 3, 3} 

GraphicalModels uses the package Graphs and a number of fundamental constructs and relationships 
associated with graphs. First we create a polynomial ring for carrying out the computations, which contains 
the entries of the covariance matrix E of a jointly normal random vector as its indctcrminatcs. Information 
about the underlying graph is stored with the polynomial ring. Hence many functions in this package take 
only ring as input, but require that it be created with gaussianRing, or markovRing in the discrete case. 

For directed acyclic graphs it is known that V^(C/ g i bai(G)) nP-D m = V(Ig) nPD m , in particular, the set 
of positive definite matrices satisfying the conditional independence constraints equals the set of covariance 
matrices in the image of the parametrization. Unfortunately, this does not imply that C/ g i b a i(G) = Ig- in 
the case of Gaussian random variables, a larger ideal, the trek ideal Tq, generated by all subdeterminants 
of the covariance matrix that vanish on the model, and satisfying C/ g i b a i(G) ^= ^G ^ Iq is sometimes equal 
to Ig (see [7]), which happens for our graph of interest G in this case. 

17 : isSubset (I , J) 

07 = true 

18 : I == J 

08 = false 

19 : J == trekIdeal(R,G) 

09 = true 

Similar computations can also be performed for graphical models with discrete random variables, and 
with other graph families. The mathematical explanation of these graphical models and their associated 
ideals appear in the remaining sections. 

2 Computing conditional independence ideals 

Conditional independence constraints on discrete or Gaussian random variables translate to rank conditions 
on certain matrices associated to the probability densities. We briefly explain these constructions here and 
how to generate these constraints in Macaulay2 using GraphicalModels. See pfl Ch. 3] for more detail. 

Let X = (Xi, . . . , X n ) be a discrete random vector where each random variable Xi has state space 
[di] = {1, 2, . . . , di}. Let d = (di, . . . , d n ). A probability distribution for X is a tensor in R dl (g) • • • §5 M. dn , 
all of whose coordinates are nonnegative and sum to one. The set of all such distributions is the probability 
simplex A^. Let pi ± ...i n = P(Xi = i\,...,X n = i n ) denote the probability of a primitive event. The 
polynomial ring in these quantities is created using the command markovRing. 

110 : d = (2,3,2); R = markovRing d 

010 : PolynomialRing 

111 : gens R 

011 = {p , p , p , p , p , p , p , p , p 

1,1,1 1,1,2 1,2,1 1,2,2 1,3,1 1,3,2 2,1,1 2,1,2 2,2,1 

P > P . P > 

2,2,2 2,3,1 2,3,2 
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For A C [n], let Xa = (x a )A be the subvector indexed by A. Let A, B, C be disjoint subsets of [n]. The 
conditional independence statement Xa -LL Xb\Xc holds if and only if the conditional distribution satisfies 

P(X A = i A ,X B = i B \X c = i c ) = P{X A = i A \X c = ic) ■ P(X B = *b\X c = i c ) 

for all iA,iB,ic- This translates into vanishing 2x2 minors of certain matrices in the probabilities p^...^. 
Those matrices are computed with the function markovMatrices, and the ideal generated by the 2x2 
minors is computed with conditionallndependenceldeal. In the following example, the two conditional 
independence statements are X\ _LL and X\ _LL X^ (:= X\ _LL Xz\X%). The ideal of vanishing minors 

has 7 quadratic generators. 

112 : s = {{{1},{2},{3», {{1},{3>,{>» 

113 : markovMatrices (R, s) 

013 = -CI p_(l,l,l) p_(i,2,i) p_(l,3,l) I, I p_(l,l,2) p_(l,2,2) p_(l,3,2) I, 

I p_(2,l,l) p_(2,2,l) p_(2,3,l) I I p_(2,l,2) p_(2,2,2) p_(2,3,2) I 
I p_(l,l,l)+P-(l,2,l)+p_(l,3,l) p_(l,l,2)+p_(l,2,2)+p_(l,3,2) 1} 
I p_(2,l,l)+P-(2,2,l)+p_(2,3,l) p_(2,l,2)+p_(2,2 ) 2)+p_(2,3,2) I 

114 : I = conditionalIndependenceIdeal(R,s) ; flatten degrees I 

014 : Ideal of R 

015 = {2, 2, 2, 2, 2, 2, 2} 

There are many different sets of conditional independence statements associated to a graph G, whose 
nodes correspond to random variables. For example, the local Markov statements of an undirected graph G 
is the set of conditional independence statements of the form X{ _LL Xy^ iuN ^y\X N ^, where N{i) is the 
set of neighbors of i in the graph G. The functions pairMarkov, localMarkov, and globalMarkov compute 
the pairwise, local, and global Markov statements, respectively, for both directed and undirected graphs. 

116 : G = graph{{l,2},{2,3},{3,4},{4,5},{l,5}}; 

117 : localMarkov G 

017 = {{{1}, {3, 4>, {5, 2», {{1, 2}, {4}, {5, 3» , {{1, 5}, {3}, {4, 2», 
{{2, 3}, {5}, {4, 1}}, {{2}, {4, 5}, {1, 3»> 

For example, the first conditional independence statement produced is X\ _LL (X3, X^)\(X2, X5). In the 
context of conditional independence, the graphical model consists of all distributions satisfying one of these 
collections of independence statements associated to the graph G. 

A Gaussian random vector, X = (Xi, . . . , X n ) ~ W(a*j £), is an rt-dimensional random vector with state 
space R n and density function 

fix) = ; 7:7- -TT7 exp [ (x — u) T Y,~ 1 (x — Lb) | , 

J{ ' (27r)"/ 2 (detE) 1 /2 y \ 2 V p; v ^' J 

where ll G M™ and E = (cr s ,t) G PD n , the cone of nxn symmetric positive definite matrices. The Gaussian 
random vector X satisfies the conditional independence statement Xa -LL X B \Xc if and only if the submatrix 
Sauc,buc := (o~ s ,t)seAuc,teBuc has rank < #C. Hence the set of all Gaussian random vectors satisfying a 
given collection of conditional independence statements naturally corresponds to a subset of PD n , and are 
studied algebraically by investigating the corresponding determinantal ideals in the polynomial ring in the 
er Sjt indeterminates. The corresponding ring is generated using the command gaussianRing. Computations 
involving conditional independence ideals with Gaussian random variables were exemplified in Section [T] 

3 Computing the vanishing ideal of a model 

The fact that graphical models can be described in two possible ways (either by a recursive factorization of 
probability distributions or by conditional independence statements) corresponds to the algebraic principle 
that varieties can be presented either parametrically or implicitly. The vanishing ideal of a model is the set 
of homogeneous polynomial relations in the joint probability distributions (for discrete random variables) or 
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in the variance- covariance parameters (for Gaussian random variables). GraphicalModels has the capability 
of computing the vanishing ideals of graphical models on directed graphs (for discrete random variables) and 
also of graphical models on directed, undirected, or mixed graphs (for Gaussian random variables). The 
vanishing ideal of an undirected graphical model for discrete random variables is a toric ideal and should be 
computed using the Macaulay2 package FourTiTwo. 

The method discreteVanishingldeal implements this capability for graphical models on discrete ran- 
dom variables. For a directed acyclic graph G on discrete random variables, the graphical model consists of 
all distributions satisfying the recursive factorization property 

p(X = i) = JJP(X W = i v \X pa{v ) = ipa(o)), 

V 

where the product runs over all vertices v of G and pa(i>) is the set of parents of v. Our implementation 
of this method does not compute the kernel of the corresponding ring map. Instead, the vanishing ideal is 
computed recursively using the factorization 

P(X = i) = P(Xi =h,.. .,X„_i = i„_i) • P(X n = i n \X pa{n) = i pa(n) ), 

where 1, . . . , n is a topological ordering of the vertices of the directed acyclic graph G. 

The following example computes the vanishing ideal of the graphical model 1— > 2 — s> 3 — > 4 on four 
binary random variables. The vanishing ideal is minimally generated by 20 quadratic binomials. 



il8 
il9 
i20 
o20 
i21 



o21 = total 




G = digraph {{1,{2}}, {2,{3»,{3,{4»,{4,{»>; 

R = markovRing (2,2,2,2); 

I = discreteVanishingldeal (R,G); 

Ideal of R 

betti mingens I 

1 

1 20 
1 . 
. 20 
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According to [3] , the vanishing ideal of a graphical model on discrete random variables is the distinguished 
component of the conditional independence ideal described by the Markov statements of the model. For the 
directed path in our previous example, the conditional independence ideal of the local Markov statements is 
a radical ideal with 3 associated primes. However, since G is a directed tree, the conditional independence 
ideal of the global Markov statements is a prime ideal and it equals the vanishing ideal of G. 



i22 
o22 
i23 



J = conditionallndependenceldeal (R, localMarkov G) ; 
Ideal of R 
I == J 



o23 = false 



i24 
o24 
i25 



K = conditionallndependenceldeal (R, globalMarkov G) ; 
Ideal of R 
I == K 



o25 = true 

The method gaussianVanishingldeal computes the vanishing ideal of a Gaussian graphical model on 
a graph, digraph, or mixed graph. If G is a mixed graph, the vanishing ideal 1q corresponds to kernel 
of the polynomial parametrization of the variance-covariance parameters o~ij in terms of the direct causal 
effect parameters Ay associated to the directed edges of G, and the noise parameters tpij associated to the 
bidirected edges of G. This parametrization is given by 

S = (7- A)- T *(7- A)" 1 , 

where £ is the variance-covariance matrix, A is the strictly upper triangular matrix with Ay = Ay if i — > j 
is a directed edge in G and otherwise, and ^ is a positive definite matrix of parameters ipij with zeros in 
each entry \l/y if there is no bidirected edge in G between i and j, and i ^ j. 
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The following example computes the vanishing ideal of the Gaussian graphical model on the mixed graph 
with directed edges {1 -J> 2,1 -» 3,2 ->■ 3,3 ->• 4} and bidirected edges {1 «-» 2,2 <-> 4}. This ideal is a 
principal ideal generated by one quartic polynomial with 8 terms. This ideal is not determinantal, i.e., it is 
not generated by the determinantal equations defining the trek ideal, which in this case is the zero ideal. 

126 : G = mixedGraph (digraph {{l,{2,3»,{2,{3}},{3,{4»},bigraph {{1,2}, {2, 4}}) ; 

127 : R = gaussianRing G; 

128 : I = gaussianVanishingldeal R; 

028 : Ideal of R 

129 : J = trekldeal (R,G) 

029 = 

An important problem in these models consists in finding which parameters are identifiable or generically 
identifiable, see [2J. The method identif yParameters can be used to solve the identihability problem 
for Gaussian graphical models on mixed graphs (also known as structural equation models). The following 
example shows that the parameter ip24 is generically identifiable by the formula ?/>24 = (013024 — cri4<723)/ci3. 

130 : H = identif yParameters R; 

131 : H#(p_(2,4))_0 

o31 =p s +s s -s s 

2,4 1,3 1,4 2,3 1,3 2,4 

In this model there are three non-generically identified parameters, identif yParameters produces a 
hash table whose entries are indexed by the parameters and contain ideals that can be used to find explicit 
rational functions for every parameter that is generically identifiable. 
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